function fsw_deep
 global N g u du
 format long e
 
 N = 1024*4; du = 2*pi/N;
 u = 2*pi*((0:N-1)/N ) ;
 k = (-N/2:N/2-1); k0 = 1*N; ki = 2; ki = 1; mu = 1e-1;
 g = 1; a = 0.2; % papa 0.28

 y(:,1) = 1 + 1.*a*exp(-ki*1i.*u);
 y(:,2) = 1.i*a*sqrt(g).*exp(-ki*1i*u);
 
 %load('N32k_tol1e10c1092951375.mat'); y = out;
 y = stokes(1.0929514,y,1e-14); 
 %c = 1.09295150 -- diverge
 %c = 1.092951375 -- converge very slowly
 
 %y(:,1) = 1 + mu*exp(-2i*u);
 %y(:,2) = -1i*mu/sqrt(5).*exp(-ki*1i*u);
 

 
 [E0 ML0 M0 AE0] = hamiltonian(y);
 cfl = 1; tmax = 10.05; 
 dt = cfl./N;
 
 t_out = 0; t_skip = 0.001; t = 0; jj = 0;
 ffilter = exp(-dt*(k/k0).^8)';
 for jjj = 1:N
  if (k(jjj) > 1) ffilter(jjj) = 0;
  end
 end
 ffilter = fftshift(ffilter);
 while t < tmax
    
    y0     = y;
    r0 = rhs(y);
    y     = y0 + 0.5*dt*r0;
    r1  = rhs(y);
    y     = y0 + 0.5*dt*r1;
    r2 = rhs(y);
    y     = y0 + dt*r2; 
    r3 = rhs(y);    
    y     = y0 + (dt/6.)*(r0  + 2.*r1  + 2.*r2  + r3 );
    
    y(:,1) = ifft(fft(y(:,1)).*ffilter);
    y(:,2) = ifft(fft(y(:,2)).*ffilter);
    
    if (t > t_out)
        figure(2)
        semilogy(k, abs(fftshift(fft(y(:,1)))), k, ifftshift(ffilter))
        axis([-N/2-10 N/2 + 10 1e-16 1e5])
        jj = jj + 1;
        [E(jj) ML(jj) M(jj) AE(jj) K(jj) P(jj)] = hamiltonian(y);
        T(jj) = t; 
        figure(4)
        semilogy(T,abs(E-E0), T, abs(M-M0), T, abs(ML) ) 
        legend('Abs. Error in energy','Abs. Error in momentum', 'Mean level', 4)
        %axis([0 tmax 1e-24 1e-2])
        figure(5)
        plot(T,K, T, P ) 
        %figure(5)
        %plot(T, K, T, P, T, K+P)
        t_out = t_out + t_skip;
        pause(0.01)
    end  
    t = t + dt;      
end
 
 
 
 
 
end

function out = Pj(in)
 global N
 k = ifftshift((-N/2:N/2-1))';
 out = 0.5*ifft((1-sign(k)).*fft(in));  
end

function out = Pj_prime(in)
 global N
 k = ifftshift((-N/2:N/2-1))'; 
 out = 0.5*ifft(1.i*(1 - sign(k)).*k.*fft(in));  
end

function out = der(in)
 global N
 k = ifftshift((-N/2:N/2-1))'; 
 out = ifft(1i*k.*fft(in));
end


function out = rhs(in)
 global g
 
 X = in(:,2).*conj(in(:,1));
 
 U = 2.*Pj(real(X(:))); 
 B = Pj_prime(in(:,2).*conj(in(:,2)));
 
 out(:,1) = 1i*(U(:).*der(in(:,1)) - der(U(:)).*in(:,1) );
 out(:,2) = 1i*(U(:).*der(in(:,2)) - in(:,1).*B(:)) + g*(in(:,1) - 1);
end

function [z_out phi_out AZ w_T Z_u] = surface_shape(in)
 global N 
 k = ifftshift((-N/2:N/2-1))';
 Z_u = find_zu(in); 
 Phi_u = -1i*in(:,2).*Z_u;
 
 z_1 = fft(Z_u)./(1i*k);  
 phi_1 = fft(Phi_u)./(1i*k); phi_1(1) = 0;
 w_T = pi*abs(k).*abs(phi_1).^2/N^2;
 
 y_1 = fft(imag(Z_u))./(1i*k*N); y_1(1) = 0;
 
 AZ = -1.*sum(abs(k).*abs(y_1).^2)*N;
 z_1(1) = 0. + 1i*AZ; z_out = ifft(z_1);
 phi_out = ifft(phi_1);
 
end

function Z = surf_p(Z_u)
 global N 

 k = ifftshift((-N/2:N/2-1))';
 z_1 = fft(Z_u)./(1i*k);
 y_1 = fft(imag(Z_u))./(1i*k); y_1(1) = 0;
 AZ = -sum(abs(k).*abs(y_1).^2)/N;
 z_1(1) = 0. + 1i*AZ; Z = ifft(z_1);

end


function [energy mean_level momentum AE Kk P] = hamiltonian(in)
 global u du g N
  [Z, Phi, AZ, Kd, z_u] = surface_shape(in);
  h = AZ/N
  %Phi_u = -1i*in(:,2).*z_u; 
  %w_T = -0.5*real(Phi).*imag(Phi_u);
  Kk = 0.5*sum(Kd);
  w_U = 0.5*g*(imag(Z)).^2.*(real(z_u));
  AE = du*sum(0.5*g*h^2.*real(z_u));
  figure(3)
  Z = surf_p(z_u);
  plot((u+real(Z)')/pi/2, imag(Z))
  axis([0 1 -0.5 0.5])
  energy = du*sum(w_U) + Kk - AE

  P = du*sum(w_U)-AE;
  mean_level = du*sum(imag(Z).*(real(1./in(:,1)))  )
  momentum   = du*sum(real(Phi).*imag(1./in(:,1)) );
end

function out = find_zu(in)
 global N
 k = (-N/2:N/2-1);
 R = in(:,1);
 Rk  = fftshift(fft(R))/N;
 rk(1:N/2+1) = Rk(N/2+1:-1:1);
 
 xk = zeros(N/2+1,1);
 xk(1)  = 1;
 for m = 2:N/2 + 1
   for n = 1:m-1
     xk(m) = xk(m) - rk(m - n + 1).*xk(n);    
   end
 end
 
 A = zeros(N/2+1,N/2+1);
 for n = 1:N/2+1
     for m = 1:n
         A(n,m) = rk(n-m+1);
     end
 end
 
 
 
 out = zeros(N,1);
 out(2:N/2+1) = xk(N/2:-1:1);
 %out(N/2 + 1) = 1;
 out = ifft(ifftshift(out))*N;

end


function out = stokes(c, in, tol)
 global g N 
 
 u = 2*pi*((0:N-1)/N )';
 k = ifftshift((-N/2:N/2-1))';
 [Z, q1, q2, q3, Z_u] = surface_shape(in);
 Z_ub = zeros(N,1);
 Z = ifft(1.i*imag(fft(Z)));
 Z_utilde = ifft( -k.*imag(fft(Z)));
 M = 0;
 while ( abs(real(M)-1) > tol )
  Z_ub = Z_utilde;  
  
  Zk = fftshift(fft(Z_utilde))/N; 
  %Zk(N/2+1) = 0.0i; 
  Z_utilde = ifft(ifftshift(Zk))*N;
  rhs = (2.*g/c^2)*Pj(imag(Z).*(Z_utilde+1));
  rhs_k = fftshift(fft( rhs ))/N;
  
  %M = sum(Zk.*Zk)./sum(rhs_k.*rhs_k) 
  M = sum(Zk.*Zk)./sum(Zk.*rhs_k)

  
  Z_utilde = rhs*M^(3/2);
  Z = surf_p(Z_utilde + 1);
  %Z_utilde = der(Z);
  b = (2.*g/c^2)
  break 
end 
 out(:,1) = 1./(Z_utilde(:) + 1);
 phi_u = c.*(Z_utilde);
 
 out(:,2) = 1i.*phi_u(:)./(Z_utilde(:)+1);
 

end



